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ABSTRACT 

The relatively large Thomson optical depth, r^-j, inferred recently from the WMAP observations suggests that 
the Universe was reionized in a more complex manner than previously believed. However, the value of t^j 
provides only an integral constraint on the history of reionization and, by itself, cannot be used to determine 
the nature of the sources responsible for this transition. Here, we show that the evolution of the ionization state 
of the intergalactic medium at high redshifts can be measured statistically using fluctuations in 21 centimeter 
radiation from neutral hydrogen. By analogy with the mathematical description of anisotropics in the cosmic 
microwave background, we develop a formalism to quantify the variations in 21 cm emission as a function of 
both frequency and angular scale. Prior to and following reionization, fluctuations in the 21 cm signal are me- 
diated by density perturbations in the distribution of matter Between these epochs, pockets of gas surrounding 
luminous objects become ionized, producing large HII regions. These "bubbles" of ionized material imprint 
features into the 21 cm power spectrum that make it possible to distinguish them from fluctuations produced by 
the density perturbations. The variation of the power spectrum with frequency can be used to infer the evolu- 
tion of this process. As has been emphasized previously by others, the absolute 21 cm signal from neutral gas 
at high redshifts is difficult to detect owing to contamination by foreground sources. However, we argue that 
this source of noise can be suppressed by comparing maps closely spaced in frequency, i.e. redshift, so that 21 
cm fluctuations from the IGM can be measured against a much brighter, but smoothly varying (in frequency) 
background. 

Subject headings: cosmology: theory - intergalactic medium - diffuse radiation 
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1. INTRODUCTION 

One of the long-standing goals of cosmology is to under- 
stand how structures have grown through time. In the usual 
paradigm, weak density perturbations were imprinted on the 
Universe during the inflationary era. These grew through 
gravitational instability, eventually forming bound halos as 
well as the cosmic web of sheets and filaments. Precise 
measurements of the cosmic microwave background (CMB) 
anisotropics have fixed the initial conditions of the picture 
(e.g., Spergel et al. 2003). The challenge now is to take struc- 
ture formation beyond the well-understood linear regime in 
order to understand how baryons collapsed into the bound 
objects that we observe today, such as galaxies and galaxy 
clusters, and to understand how these objects affect their sur- 
roundings. At low or moderate redshifts (z < 6), galaxies and 
quasars can be studied in detail with existing technology. Un- 
fortunately, the first generations of protogalaxies are not yet 
accessible observationally. Their properties are nevertheless 
crucial to understanding both later generations of galaxies, 
which form out of these early protogalaxies in any hierarchi- 
cal picture of structure formation, and the gross evolution of 
baryons in the Universe, because these objects exhibit strong 
feedback on their surroundings. Perhaps the most important 
such channel is the reionization of the intergalactic medium 
(IGM). When the first protogalaxies or quasars form, they ion- 
ize pockets of surrounding gas. These H II regions grow with 
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time and eventually overlap. The timing, morphology, and 
duration of this event contain a wealth of information about 
both the ionizing sources and the IGM (e.g., Wyithe & Loeb 
2003; Cen 2003; Haiman & Holder 2003; Mackey et al. 2003; 
Yoshida et al. 2003a,d). 

A great deal of effort has gone into constraining the tran- 
sition from a neutral to ionized IGM. Unfortunately, exist- 
ing observational techniques are not optimized to the needed 
measurements; they have provided tantalizing constraints on 
reionization but cannot be used to map the event in detail. The 
most straightforward method is to extend the "Lya forest" to 
high redshifts: regions with relatively large H I densities ap- 
pear as absorption troughs in quasar spectra, which presum- 
ably deepen and come to dominate the spectra as we approach 
the reionization epoch. Indeed, spectra of z ~ 6 quasars se- 
lected from the Sloan Digital Sky Survey"^ (SDSS) show at 
least one extended region of zero transmission (Becker et al. 
2001), indicating that the ionizing background is rising at this 
time (Fan et al. 2002). However, the optical depth of the 
IGM to Lya absorption is Tiya « 6.45 x 10%[(1 +z)/10]^/2 
(Gunn & Peterson 1965), where xh is the neutral fraction. A 
neutral fraction xh > 10^^ will therefore render the absorp- 
tion trough completely black; quasar absorption spectra can 
clearly probe only the latest stages of reionization. 

A second constraint comes from the effects of the ionized 
gas on the CMB. The free electrons Thomson scatter the CMB 
photons, washing out the intrinsic anisotropics but generating 
a polarization signal. The total scattering optical depth t^s is 
proportional to the column density of ionized hydrogen, so it 
provides an integral constraint on the reionization history. Re- 
cently, the Wilkinson Microwave Anisotropy Probe^ (WMAP) 
used the polarization signal to measure a large Tes, indicating 

■* See http://www.sdss.org/. 
' See http://map.gsfc.nasa.gov/. 
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that reionization began at Zr ?J 14 (Kogut et al. 2003; Spergel 
et al. 2003). More detailed information on the reionization 
history could be obtained by measuring the (large) angular 
scales over which CMB polarization is generated (Zaldar- 
riaga 1997; Kaplinghat et al. 2003; Hu & Holder 2003) or the 
(small) scales over which secondary anisotropics are gener- 
ated by the patchiness of reionization (Gruzinov & Hu 1998; 
Knox et al. 1998), but these signals promise to be difficult to 
extract (Holder et al. 2003; Santos et al. 2003). A third con- 
straint comes from measurements of the temperature of the 
Lya forest at z ~ 2-4, which suggest an order unity change 
in the ionized fraction at Zr ^ 10 (Theuns et al. 2002; Hui & 
Haiman 2003), although this argument depends on the timing 
and history of He II reionization (e.g., Sokasian et al. 2002). 

Taken together, these three sets of observations imply a 
complex reionization history extending over a large redshift 
interval (Az ^ 10). This is inconsistent with a "generic" pic- 
ture of fast reionization (e.g., Barkana & Loeb 2001, and ref- 
erences therein). The results seem to indicate strong evolution 
in the sources responsible for reionization, and a detailed mea- 
surement of the reionization history would contain a rich set 
of information about early structure formation (Sokasian et al. 
2003a; Wyithe & Loeb 2003; Cen 2003; Haiman & Holder 
2003). The optimal reionization experiment would: (1) be 
sensitive to order unity changes in xh (to probe the crucial 
middle stages of reionization), (2) provide measurements that 
are well-localized along the line of sight (rather than a single 
integral constraint), and (3) not require the presence of bright 
background sources, which may be rare at high redshifts. The 
most promising candidate proposed to date is the 21 cm hy- 
perfine transition of neutral hydrogen in the IGM (Field 1958, 
1959a), which fulfills all three of these criteria. So long as the 
excitation temperature 7^ of the 21 cm transition in a region 
of the IGM differs from the CMB temperature, that region 
will appear in either emission (if Ts > Tqmb) or absorption (if 
Ts < Tcmb) when viewed against the CMB. Variations in the 
density of neutral hydrogen (due either to large-scale struc- 
ture or to H II regions) would appear as fluctuations in the sky 
brightness of this transition. Because it is a line transition, 
the fluctuations can also be well-localized in redshift space. 
Thus, in principle, high resolution observations of the 21 cm 
transition in both frequency and angle can provide a three- 
dimensional map of reionization. Together with radio absorp- 
tion spectra of bright background sources (which can probe 
much smaller physical scales in the IGM; CarilU et al. 2002; 
Furlanetto & Loeb 2002), these observations promise to shed 
light both on the early growth of structure and on reionization. 

The physics of this transition has been well-studied in the 
cosmological context. Early work focused on fluctuations 
due to large-scale structure (Scott & Rees 1990; Kumar et al. 
1995; Madau et al. 1997; Tozzi et al. 2000; lUev et al. 2002), 
because the signals could be estimated through linear cosmo- 
logical perturbation theory. Shaver et al. (1999) were the first 
to explicitly consider the signal at reionization, although they 
focused on the "all-sky" signature rather than the fluctuations. 
Recently, Ciardi & Madau (2003) and Furlanetto et al. (2003) 
used numerical simulations of reionization to estimate how 
the fluctuations would behave during that epoch. We show 
in Figure 1 three time slices from the simulation analysis de- 
scribed by Furlanetto et al. (2003), corresponding to the early, 
middle, and late stages of reionization (from left to right). It 
is clear that both the mean signal and the fluctuations drop 
abruptly. Interestingly, the fluctuations during reionization 
have a very different morphology than those due to large-scale 



structure; the spectrum of fluctuations thus has the potential to 
constrain the process of reionization. 

In this paper, we present a new approach to 21 cm fluc- 
tuations. We draw an analogy between these measurements 
and those of the CMB: in both cases we wish to measure the 
level of inhomogeneity as a function of scale on the sky. Pre- 
vious treatments of the 21 cm signal have focused on mea- 
suring fluctuations on a particular patch of the sky, implicitly 
referring to imaging observations. Here we show that a sta- 
tistical treatment of the fluctuation power spectrum contains a 
great deal of information about reionization. Furthermore, a 
large set of tools for CMB predictions and data analysis has 
already been developed (see Hu & Dodelson 2002 for a re- 
view), so there is much to be gained by connecting the two. 
Indeed, some steps in this direction have already been taken 
by Pen (2003) and Cooray (2003), both of whom considered 
the effects of lensing on the 21 cm signal. The analogy with 
the CMB is not perfect, however, because the 2 1 cm signal 
can be separated in redshift space; in other words, we can 
make maps at a series of frequencies, each of which samples 
an independent volume. In this sense the analogy is closer to 
redshift surveys (e.g., Peebles 1980) or weak lensing tomog- 
raphy (e.g., Hu 1999). We therefore develop our formalism 
with explicit consideration of how multifrequency informa- 
tion can be used with power spectrum statistics, in effect gen- 
eralizing the methodology used to analyze CMB anisotropics. 
After reviewing the physics of the 21 cm transition in §2, we 
show how to compute the angular power spectrum of 21 cm 
fluctuations in §3. 

In §4 we show some simple apphcations of our approach. 
We give predictions for the angular power spectrum of 21 cm 
fluctuations from a fully neutral medium and for a simple toy 
model of reionization. In the former case, fluctuations in the 
signal are due orily to large-scale structure. We show that in 
this regime 21 cm measurements essentially yield the power 
spectrum of density fluctuations (see also Pen 2003). We then 
show that variations in the neutral fraction during reionization 
distort the power spectrum. 

Another advantage of our approach is that the angular 
power spectra are closely related to the physically observed 
quantities. This is especially true for the interferometers that 
will most likely be used to measure the redshifted 21 cm sig- 
nal. Our results thus connect theoretical predictions to the 
potential observations. For example, inhomogeneities in the 
21 cm signal must be separated from fluctuations in any fore- 
ground sources. This is particularly important because the 
absolute foreground signal will swamp the 21 cm signal by 
many orders of magnitude. While Galactic foregrounds are 
expected to be fairly smooth on the relevant angular scales, 
faint radio galaxies, starbursts, and even the galaxies respon- 
sible for reionization fluctuate strongly on arcminute scales 
and dominate those of the 21 cm signal by at least an or- 
der of magnitude (Di Matteo et al. 2002; Oh & Mack 2003). 
These results have been used to argue that the prospects for 
large angular-scale measurements of the reionization epoch 
are dim. However, both Di Matteo et al. (2002) and Oh & 
Mack (2003) also pointed out that all of the (known) fore- 
ground sources have featureless power-law spectra. Both sug- 
gested that the foregrounds could therefore be removed in 
frequency space. As an example, consider the simple case 
in which every foreground source has the same spectral in- 
dex. Then the foregrounds between two maps at nearby fre- 
quencies would be exactly correlated, while the 21 cm fluctu- 
ations wiU be uncorrelated because each frequency samples 
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Fig. 1 . — The brightness temperature of the 21 cm transition at several redshifts, as predicted by the "late reionization" simulation analyzed in Furlanetto et al. 
(2003). Each panel corresponds to the same shce of the simulation box (with width lO/j^' comoving Mpc and depth Aj/ = 0.1 MHz), at z = 12.1, 9.2, and 7.6, 
from left to light. The three epochs shown correspond to the early, middle, and late stages of reionization in this simulation. (For details about the simulations, 
see Sokasian et al. 2001; Springel & Hemquist 2003a,b.) 



an independent volume of the IGM. Comparing two maps 
closely spaced in redshift therefore allows one to remove the 
foreground component. We show in §5 how the foreground 
sources can be modeled with our multi-frequency formalism. 
In §6 we quantify how well their contamination can be re- 
moved. We find that foregrounds are much less important than 
previously assumed so long as the range of allowed spectral 
indices for faint sources is similar to that already measured 
for brighter sources. Finally, we estimate in §7 how well the 
power spectrum can be measured with the next generation of 
low-frequency radio telescopes, and we conclude in §8. 

When necessary, we assume a ACDM cosmology with 
fl„, = 0.3, Q-A = 0.7, rtb = 0.04, Ho = lOO/z km s"' Mpc"' 
(with h — 0.7), and a scale-invariant primordial power spec- 
trum with n = 1 normalized to ag = 0.85 at the present day. 

2. 21 CM RADIATION FROM THE INTERGALACTIC MEDIUM 

The optical depth of a patch of the IGM in the hyperfine 
transition is (Field 1959a) 
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Here vq — 1420.4 MHz is the rest-frame hyperfine transition 
frequency, Ajo = 2.85 x 10^''' s^' is the spontaneous emis- 
sion coefficient for the transition, is the spin temperature 
of the IGM (i.e., the excitation temperature of the hyperfine 
transition), Tcmb = 2.73(1 +z) K is the CMB temperature at 
redshift z, and nni is the local neutral hydrogen density. In the 
second equality, we have assumed sufficiently high redshifts 

such that H{z) ~ Hofll^^{l +z)^^^ (which is well-satisfied 
in the era we study, z > 6). The local baryon overdensity 
is 1 + (5 = p/p and xh is the neutral fraction. The radiative 



transfer equation in the Rayleigh- Jeans limit then tells us that 
the brightness temperature of a patch of the sky (in its rest 
frame) is Th = Tcmb^^^ + 75(1 — e^'^). We define 6T{v) to be 
the observed brightness temperature increment between this 
patch, at an observed frequency v corresponding to a redshift 
1 + z = vq/v, and the CMB: 
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Assuming that the radiation background includes only the 
CMB, the H I spin temperature is (Field 1958) 



TcMB + + jLya ^L; 
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The second term describes collisional excitation of the hyper- 
fine transition, which couples Ts to the gas kinetic temperature 
T]^. The coupling coefficient is 



_ Cio 7* 



(4) 



where Cio(7a:) oc hh is the collisional de-excitation rate of the 
(higher-energy) triplet hyperfine level (Allison & Dalgarno 
1969) and 7; = Inhiyo/k = 0.068 K. For 7> - 1000 K, the 
coupling becomes strong when 1+5>5[(1+ z)/20]^. The 
third term in equation (3) describes the Wouthuysen-Field ef- 
fect, in which Lya pumping couples the spin temperature to 
the color temperature of the radiation field Ti,ya (Wouthuysen 
1952; Field 1958). We note that TLy^ = Tk so long as the 
medium is optically thick to Lya photons (Field 1959b). Es- 
sentially, the dipole selection rules allow a transition between 
the two hyperfine levels of the ground state mediated by the 
absorption and subsequent re-emission of a Lya photon. The 
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excitation and de-excitation rates are then controlled by the 
color temperature of the radiation field near the line center, 
which (for a sufficiently large number of scatterings) must be 
in thermodynamic equilibrium with the gas temperature. The 
coupling constant for this process is 



while a given angular scale A9 corresponds to 



Pio n 



(5) 



^10 Ti^ya 

where Piq is the indirect de-excitation rate of the triplet level 
due to absorption of a Lya photon followed by decay to the 
singlet level. For a diffuse Lya background, Madau et al. 
(1997) showed that 

Pio « 1.3 X 10-'V_2i s-S (6) 

where 7_2i is the intensity of the background radiation field at 
the Lya frequency in units of 10"^' erg cm^^ s~' Hz~' sr~'. 
Lya pumping effectively couples Ts and Tk when 7_2i > 1. 
Ciardi & Madau (2003) argue that J-21 > 1 even at z > 20. 
If so, Ts ~ Tk throughout the diffuse IGM, even though the 
densities are well below the threshold for collisional couphng. 

The spin temperature and optical depth therefore depend 
on the kinetic temperature of the IGM. Once Thomson scat- 
tering of CMB photons becomes inefficient at the thermal de- 
coupling redshift Zd ~ 140, the IGM cools adiabatically until 
the first objects collapse (Couchman & Rees 1986). During 
this era, Tk < Tqmb- The cooling trend reverses itself as soon 
as significant structure begins to form, but the subsequent 
temperature evolution is both inhomogeneous and highly un- 
certain. While early estimates suggested that Lya photons 
themselves would inject significant thermal energy into the 
IGM, Chen & Miralda-Escude (2003) showed that this heat- 
ing channel is in reality quite slow. Instead, X-rays (primar- 
ily from supernovae or accreting black holes) and shocks are 
Ukely to control the temperature evolution of the IGM. We 
expect shocks to heat overdense structures like sheets, fila- 
ments, and virialized halos to Tk > Tcmb and radiative feed- 
back from stars and quasars to heat the rest of the gas. Most 
estimates suggest that the two processes will rapidly heat the 
IGM to Tk > Tcmb (Venkatesan et al. 2001; Chen & Miralda- 
Escude 2003; Gnedin & Shaver 2003). The topology of the 
two cases of course differs; shock heating will tend to exag- 
gerate brightness temperature differences by separating warm, 
dense regions from cool voids, while X-ray heating will in- 
duce a smooth temperature distribution. 

In developing our formalism, we allow 5T{i/) to depend on 
the parameters S,xh, and T^. This is the most general case 
possible (once the cosmological parameters are fixed) and al- 
lows one to incorporate the full range of physics when neces- 
sary. However, the arguments above suggest that the situation 
most relevant to observations has Ts ^ Tk ^ Tcmb- In this 
limit, the temperature factor in equation (2) approaches unity 
and the signal is independent of Ts- Thus we essentially as- 
sume an era of significant X-ray heating; in this scenario, the 
extra heating in shocked dense regions can be ignored. For 
simplicity, we will restrict ourselves to this case for the illus- 
trative examples in §4. We emphasize that this is, however, 
an important assumption. If, for example, significant heating 
does not occur until the early stages of reionization, STi, will 
have a much more complicated distribution than we consider 
here. 

Finally, to orient the reader, we note that an observed band- 
width Au corresponds to a comoving distance 

L«1.7(,r.=^l(^|'^75£^ '^'mpc (7) 
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over the relevant redshift range. 

3. BASIC FORMALISM 

We now show how to compute the angular power spectrum 
of 21 cm fluctuations. Unfortunately, a given patch of the sky 
observed with frequency bandwidth Av does not correspond 
directly to a physical volume of the Universe because the ob- 
servation is performed in redshift space: peculiar velocities 
can move a parcel of gas into or out of this channel. However, 
redshift space distortions will be unimportant if Av jv > v/c 
where v is the typical random bulk velocity of the gas. The 
importance of redshift space distortions is determined by the 
ratio 
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For large values of TZ, redshift distortions have only a 
marginal effect. Furlanetto et al. (2003) have shown that, 
for typical survey geometries, redshift space distortions am- 
phfy the signal by at most ~ 25% (see also Tozzi et al. 2000). 
For simplicity we will ignore them in what follows. In future 
work, we will examine the significance of redshift distortions 
in detail using numerical simulations and asses what cosmo- 
logical information can be extracted from their detection.^ 

If redshift distortions can be neglected, there is a one-to-one 
correspondence between frequency and redshift. The band- 
width of the experiment is characterized by some response 
function W{y). The observed brightness temperature is of the 
form 

T{n,y) = Toiro) J W„(r)V(n,r), (10) 

where n is the direction of observation, Tq is a normaliza- 
tion constant which depends on redshift and ^'(n, r) is the di- 
mensionless brightness temperature, ip{n,r) = {1 + 6)xh{Ts — 
Tcmb)/Ts. Note that ST{y) = To{ro)'ilj{n,r) in equation (2). 
The projection window Wro{r) is a function peaked at ro, the 
radial distance corresponding to the observed frequency u, 
and has a width Sr. 
We can expand ^p{x) as a Fourier series. 



V(x) = 
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where ji{x) are the spherical Bessel functions and Yi,„{h) are 
the spherial harmonics. The statistics of V' are determined by 
its power spectrum,^ 

(V'(ki)^(k2)) = (2^)'(5^(ki +k2)P^iki). (12) 

* In principle, if one knows the dark matter power spectrum, comparison 
with the observed spectrum could allow one to extract the peculiar velocities 
directly. This procedure would be easiest when fluctuations in xh can be 
ignored, i.e. in the very early stages of reionization. 

' Actually, ip is unlikely to be a true Gaussian random field because of the 
distribution of xh, so its statistics are not entirely determined by (k). 



5 



(13) 



We use equations (10) and (11) to calculate the spherical 
harmonic decomposition aim of the observed temperature: 

■i'f- 

J (271 

ai{k,u)=To{ro) J drWr,{r)ji{kr). 
We can define the angular power spectrum as 

j^-^P^{k)ai{k,vi)ai{k,V2) 

/dk 
—Al{k)ai{k,ur)ai{k,U2)iU) 

Here {k) = k^P^ / (27r^) and we have used the isotropy of 
P^. This formula encodes both the case of the power spectrum 
of maps at one particular frequency (when v\ = 1^2) as well as 
the correlations between maps at different frequencies. Equa- 
tion (14) forms the basis for the analysis that makes it possible 
to relate the 21 cm fluctuations to the evolution of the ioniza- 
tion state of the IGM. Note that the C/'s approach zero as 
and 1^2 depart from each other for two reasons: because the 
frequency-space window functions no longer overlap in this 
Umit and because the fluctuations in the IGM are uncorrelated 
on large scales. The latter property can be used to separate the 
21 cm signal from contamination by foreground sources that 
vary smoothly with frequency. 

We can investigate the behavior of equation (14) by con- 
sidering various limits. We want to understand how equa- 
tion (14) depends on the width of the response function, 5r 
(which describes the bandwidth of the observation). The fig- 
ure of merit is I6r/r (i.e. the ratio of the radial to transverse 
scales probed by the observation). We first consider the limit 
in which the response function can be considered to be a delta 
function, Idr/r <^ 1. In that case. 
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so that 
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In the limit in which the power spectrum can be approximated 
by a power law, A^(A;) = {k/k^)" we then have 
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With this definition f{n,l) is a very weak function of both / 
and n and it is normahzed so that f{Q,l) = 1. Thus in this 
Umit {I6r/r <^ 1), the temperature fluctuations simply trace 
the underlying tp fluctuations. 
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can approximate the integral in a different way, 
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We now consider the opposite regime, I5r/r^ 1, the large 
bandwidth or Limber Umit (Limber 1953; Peebles 1980). We 



The Bessel functions ji (x) are very small for x <l and start 
to oscillate when I. Thus, the integral over k will receive 
contributions only from a region around k^ljr with width of 
order Ak ~ Xjbr. Modes with k>k-\- Ak will be out of phase 
for typical separations of the two points r\ and ^2. In this 
regime we can approximate P^(K) k, P^{1 j r\). The integral 
over k is then proportional to b^(r\ — ¥2) /r\ so that 

Q{yuy2)=T^ j drW\r)^^ 

a T, ^^{llro)j^^. (20) 

This is the standard Limber's equation, widely used in the 
context of weak lensing (e.g. Kaiser 1992). 

Equations (18) and (20) are easy to understand. For a suf- 
ficiently narrow frequency response {I6r/r <C 1), the angu- 
lar fluctuations directly trace those of the underlying if} field. 
However, for a surface of finite width 5r only those modes 
with radial k<l/Sr can contribute because the response func- 
tion averages out larger k modes. Thus when the surface be- 
comes too thick in the radial direction, angular fluctuations 
are no longer of order ~ k^P{k) but become ^ k^P{k)/Sr, as 
seen in equation (20). 

To estimate the / at which the width of the surface begins to 
damp the fluctuations, we can consider an Einstein de-Sitter 
universe in which a(T) = (t/tq)^, where a is the scale factor 
and r is the conformal time. In this case 6t/t = {\/2)Sa/a = 
(1/2) Ai// u. Thus, the changes in radial distance r = tq — r 
are just 

,dr, 1 Av 

where ly is the observed frequency, i' = 1^0/ {1 + z). For Ai' = 
0.2 MHz the corresponding value of / from equation (21) is 
I ^ 5000, or arcminute scales. The damping of fluctuations 
for large bandwidths has important impUcations for the choice 
of ( Ai/, /) in a given observation (see Figure 4 below and the 
discussion thereof). 

4. SIMPLE MODEL FOR THE CORRELATIONS OF THE 
BRIGHTNESS TEMPERATURE 

In this section we will make a simple model for the correla- 
tions of the dimensionless brightness temperature, ^. A more 
detailed study using simulations is left for future work. 

We assume that Ts > Tcmb so that V = :ch(1 + To cal- 
culate the power spectrum of tp we need a model for the cor- 
relations of the neutral fraction, xh- We will follow a similar 
treatment as that used to model the effect of patchy reioniza- 
tion on CMB anisotropics (e.g. Gruzinov & Hu 1998; Knox 
et al. 1998). For simplicity, we wiU model the the fluctuations 
in Xh as if they were produced by a set of uncorrelated "bub- 
bles" of typical size R. We denote the average value of xh as 
Xh- We allow both R and xh to depend on redshift. Under this 
simpUfying assumption we can model the correlations as 

{Xh{Xi)Xh{X2)) = 4 + {xh - 4)/(^i2/^), (22) 
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where f{x) is a function with the following limits: f{x) w 1 
for X <C 1 and f{x) k, for x » 1. The details of this func- 
tion depend on the shape of the bubbles. If the bubbles are 
correlated, one should interpret R as an effective size that de- 
pends on the correlation length between distinct bubbles. To 
calculate observables, we will take f{x) — exp(— 

The correlations of i/', i^i^n) = (?/'(xi)?/'(x2)) — (V')^^ be- 
come 

^i{xn) = [xh + {xH-XH)f{xn/R)]£.{xn) + 

{xH~x]j)f{xn/R) + ?7(xi2)[2xh + ri{xn)] , (23) 

where ^{xn) — ((5(xi)(5(x2)) is the correlation function of the 
density field and rj{x\2) — ((5(xi)x/f(x2)) gives the cross cor- 
relation between the density and neutral fraction fields. To 
keep things as simple as possible we will ignore this last term 
in what follows. We will explore the consequences of includ- 
ing it in future work. The correlation function is related to the 
density power spectrum by 



(24) 



sinfer 
IT' 

Equation (23) has the following hmits: 

tx{xi2)r^XH£,{xi2)-¥{xH-x]f) {xn<.R), 
tJ-ixi2)~xl^ixi2) (xi2>-R). (25) 

When xi2 ^ R both points are basically independent as 
far as the correlations in xh are concerned, so the correla- 
tions are given by those of the density field times the prob- 
ability that each of the two points falls in a neutral region, 
xjj. On scales smaller than the "bubble" size, both points 
fall either inside or outside a "bubble" so only one factor of 
Xh multiplies ^. On top of the fluctuations produced by the 
density there are those created by the presence of the bub- 
bles, (xh — xjj). To illustrate the behavior of equation (23) 
we can take the correlation function of the density to be a 
power law, ^(x) = (x/xq)"". Moreover, assume that xq < R 
and Xh ~ 0.5. On scales smaller than xo, /i(x) « xh(x/xo)^". 
In the range xq < x < R, /i(x) « xh{1 —xh), while forx ^ R, 
/i(x) « x^(x/xo)^". Thus there is a feature in the correlation 
function on the scale of the bubbles. In this simple model the 

correlations trace the matter correlations on both large and 
small scales (but with different amplitudes). On scales inter- 
mediate between the size of the bubbles and the non-linear 
scale the correlation function flattens out. 

We can Fourier transform equation (23) (with 77 = 0) to ob- 
tain an expression for the power spectrum of ip, 

+ ixH-xl)A\{k), 
where we have introduced 



(26) 



AMk) = 



k\f{k) 
2^2 ' 

2^2 



(2^ 



Pp(k-k')/(k'), (27) 



with f{k) the Fourier transform of /(x) and ^'p(k) is the power 
spectrum of the density fluctuations. 

In Figure 2 we show the power spectrum of ij: calculated 
from equation (26). For illustrative purposes, we used a sim- 
ple model for the mean neutral fraction as a function of red- 
shift, 

l+exp[-(z-Zo)/Az] 




k[hMpc-'] 



Fig. 2. — Power spectrum of ip fluctuations, = k^P^ jl-K-, in tlie model 
described in tlie text. The lines coirespond to z = 12, 11, 10,9, and 8. At 
z = 10 (solid black line), xh = 1/2. The long-dashed blue lines con'espond 
to z = 12 and 11 at the beginning of reionization and the short-dashed red 
lines to z = 9 and 8. The dotted green line shows the power spectrum of 
the density field at z = 10. Note that we use the linear dark matter power 
spectrum in computing the density fluctuations. 



with zo = 10 and Az = 0.5. For the correlation length R 
we used a maximum value of 3/!~'Mpc when xh = 0.5 and 
smaller values when xh deviates in both directions from 0.5 
so as to keep the number density of bubbles fixed. We 
show results for, {z,xh,R) = (12,0.98,1.24), (11,0.88,2.24), 
(10,0.5,3), (9,0.12,2.24), (8,0.02,1.24). For comparison, we 
also show the power spectrum of matter fluctuations at red- 
shift z = 10. The form of this choice for the evolution of 
the mean neutral fraction is motivated by numerical simula- 
tions of reionization (e.g. Figures 5 and 9 of Sokasian et al. 
2003a,b, respectively). Note that we have simply used the lin- 
ear dark matter power spectrum at the appropriate redshift in 
computing for the Figure. In reality, of course, we should 
use the full gas power spectrum, which includes the nonlinear 
growth of structure on small scale and smoothing due to the 
finite pressure of the gas. The differences are, however, not 
large (see, for example. Figure 2 of Furlanetto et al. 2003), so 
the linear dark matter spectrum will suffice for our purposes 
here. We will consider modifications due to the true power 
spectrum in future work. 

Before reionization begins, the power spectrum of ip is sim- 
ply the power spectrum of the matter fluctuations. As the 
neutral fraction decreases it develops a feature on the scale 
of the bubbles, roughly atk^2/R. The feature has maximum 
amplitude when xh = 0.5. As the neutral fraction decreases 
further, the 21 cm fluctuations disappear because there is no 
longer any more neutral hydrogen to emit radiation. Note that 
when the bubbles appear, the power spectrum on the largest 
scales is a power law A^, cx k^, corresponding to Poisson fluc- 
tuations. That is, the Poisson fluctuations induced by the dis- 
crete nature of the bubbles dominates over the fluctuations 
resulting from 5. 

In Figure 3 we show the corresponding angular power spec- 
trum ST = [l{l + l)C//27r]'/2 with a window function W,„(r) 
of Gaussian shape, centered at the appropriate redshift and 
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Fig. 3. — Angular power spectrum of 21cm fluctuations at z = 12, 1 1, 10, 9, 
and 8 in the model described in the text (the curves are the same as in Figure 
2). The bandwidth is 0.2 MHz. 



with a full width at half maximum of 0.2 MHz. The angular 
power spectrum traces the behavior of the 3D power spectrum 
shown in Figure 2. It develops a feature on the scale of the 
bubbles, I ^kr ^ 2r/R. We emphasize that the precise shapes 
of the curves like those in Figures 2 and 3 depend on the mor- 
phology of the ionized regions. Here, we have adopted a sim- 
ple model in which these "bubbles" are described by spheres 
whose size evolves simply with redshift. In reality, the ion- 
ized regions have complicated morphology and evolution, as 
indicated by Figure 1, that will likely imprint more complex 
features into the power spectra than suggested by Figures 2 
and 3. We will examine these issues further in due course. 

On small scales, the angular power spectrum in Figure 3 
changes slope and stops tracing A^. This occurs when the 
window in frequency becomes too wide and we enter into the 
Limber regime [equation (20)], where the angular fluctuations 
are damped by one power of k. To illustrate this further, we 
plot the angular power spectrum for several spectral widths 
in Figure 4. As the width of the filter increases, the level of 
fluctuations decreases. The Figure shows that the damping 
is insignificant until the Limber regime is reached at I5r/r ~ 
1 ; in the regime where the fluctuations are dominated by the 
bubbles, this happens when 5r/R ^ \. For bandwidths larger 
than this, the level of fluctuations scales as l/(5r c>c 1/ for a 
fixed angular scale. On the other hand, the signal per channel 
is proportional to the bandwidth. Thus, if one is interested in 
a particular angular scale Z, it is best to choose the bandwidth 
Av such that I5r/r < I. 

5 . MODELING THE CONTAMINANTS 

One of the major challenges in observing this signal is sep- 
arating it from the many other (stronger) low-frequency ra- 
dio sources. There are several potential sources of noise, 
but those associated with foreground sources are likely to be 
smooth in frequency space. In this section we will show how 
this smoothness can be used to remove this contamination. 
We will explicitly consider point source foregrounds, as dis- 
cussed by Di Matteo et al. (2002) and Oh & Mack (2003), 
but our method can be easily be extended to other types (such 




10 100 1000 



Fig. 4. — Top panel: Angular power spectrum at redshift z = 10 for obser- 
vations done with filters of full width at half maximum (0.1,0.2,0.4, 0.8, 1.6, 
3.2, 6.4) MHz. The dotted line in the top panel shows the power spectrum of 
(5 at z = 10. Bottom panel: Ratio between the power spectra in the top panel 
and one corresponding to a delta function filter. 



as diffuse Galactic synchrotron emission), provided that their 
power spectra can be estimated. 

Let us assume that there is a collection of different types of 
sources with a range of spectral indices C- We describe them 
by a luminosity function d^n/dSdC, which gives the average 
number of sources per steradian per unit flux S and spectral 
index C,. We will assume that the sources are clustered but 
not necessarily that sources of different spectral indices are 
perfectly correlated. The clustering on the sky is described by 
the angular power spectrum. 

Km (Ci Km (C2)) = S„J„nmq, (Cl , C2) (29) 

= s,,i,6„„„,ci (Cl , C2)y^q(Ci)q(C2), 

where C[(Ci,C2) is the cross correlation between the pop- 
ulations and q(C) = C[(C,C) is the auto correlation of a 
given population. These are nothing more than the Legen- 
dre transforms of the respective correlation functions. We 
have also introduced the correlation coefficient Cf {(1,(2) = 

q(Ci,C2)/Vq(Ci)q(C2). 

The power spectrum of fluctuations on the sky produced by 
these sources is given by, 

c,(.i.2)^/...C^.^(^)^(f)^ ^ ^-) 
xC,(Ci,C2)y^y^(^j [j) . 

The first term is the Poisson contribution and the second 
comes from the clustering of_the sources. We have denoted 
the average spectral index by 

We will eventually show that what interferes with a mea- 
surement of the 21 cm signal is the fact that the foreground 
maps in different frequencies may not be perfectly correlated. 
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If the maps were perfectly correlated, then one could in some 
sense subtract the map at one frequency from another and 
clean the map of all contamination. In our simple model the 
presence of sources with several spectral indices is at the heart 
of the fact that maps at different frequencies become uncorre- 
lated. 

To make some progress, we adopt some assumptions about 
the different quantities that enter into equation (30). We will 

'"i/ic)^^''";^;'^ (31) 



dSdC dS-^ dS V2^SC 

where S( measures the range of spectral indices of the 
sources. This Gaussian form is a reasonable description of 
the spectral index distribution of low-frequency radio sources 
measured by Cohen et al. (2003). We need to model the cor- 
relation coefficient between different sources. We know that 
(C> C) = 1 and should decay as the two ('& depart from each 
other. We will take 

where ct^ measures how sources with different spectral in- 
dices become uncorrelated. For simplicity, we will assume 
that C is independent of /. This could be relaxed easily but 
would make our expression a bit more complicated. To the 
extent that all sources are tracing the same underlying distri- 
bution of matter, they should be perfectly correlated, even if 
their bias is somewhat different. Only the stochastic part of 
the bias contributes to the loss of correlation. Thus we expect 
£7^ to be large compared to SC that is, all the sources should 
be well correlated on the sky. Finally, we will approximate 



/q(C) 



1 



(33) 



yq(C) 2flflnC^C 

meaning that we will compute quantities of interest only to 
lowest order in the change of clustering with population. 

We use the above formulae to calculate the correlation co- 
efficient between maps at different frequencies, 



Il{v\,l'2) = 



(34) 



\JCi{vi,vi)Ci{v2,V2) 
We calculate //(j^i,^'2) as a series in \n{v\/i'2) and to lowest 
order in JlnQ/Jln^and^C/o'C and obtain 



1-1-7 (1+7)^ Vo-f 



d\nC 



7 



1+27(1+7) 



4(1+7)2 



(35) 



2(1+7)3 

where 7 measures the importance of the Poisson term relative 
to the term due to clustering: 



7 = 



^poisson 

^ 

'duster ' 



Cr"™"= / dS—S^, 



-^cluster 
-I 



dS 

r2 



These are the standard formulae for the Poisson and clustering 
contributions for a single population of sources (e.g. Peebles 
1980). Equivalently from equation (30) and to lowest order in 
SC/(7(^ andJlnQ/Jln^wefind 

Ci{i?,u)=Cf°'''"" + q{C)I^. (37) 

We can point out a few interesting things about equation 
(35). First of all, the departure from unity is proportional to 
S(^ln^{iyi/i^2)- As expected, there is a loss of correlation only 
to the extent that there are sources with a variety of spectral 
indices in the mix. There are several contributions to the loss 
of correlation. The first term on each line (proportional to 
7/[l + 7]) is a direct consequence of the Poisson part; they 
go to zero as 7 ^ 0. They occur because, if the Poisson con- 
tribution is dominant, there is a chance of getting a different 
spectral index in different regions of the sky. That is to say, 
there are a few sources in any given mode on the sky and 
so the fluctuations in the spectral index of actual sources that 
happen to be in each region will lead to patterns on the sky 
that are not identical at different frequencies. The other terms 
come from the clustering part (they go to zero as 7 — > 00). 
Those terms only appear if different sources are not perfectly 
correlated (terms proportional to S(/(j^) or if they cluster dif- 
ferently (terms proportional to d\nCi/d\n(). 

The main point of equation (35) is that the difference be- 
tween the cross-correlation and unity scales as In^ {vi jv2), 
so it should be quite small. This will imply that the noise from 
the smooth foregrounds is unimportant. 

6. THE EFFECT OF CONTAMINATION 

In this section we will show how the frequency information 
can be used to discriminate the 21 cm signal from sources of 
contamination. The basic point is that while the contaminants 
are presumed to be smooth as a function of frequency, the 2 1 
cm signal varies very rapidly. A small change in frequency 
of order a fraction of a MHz is already enough to sample an 
effectively different part of the Universe. In what follows, we 
will present two derivations that show that the measurement 
is only contaminated by the parts of the foregrounds that are 
uncorrelated between neighboring frequencies. For simplicity 
we will only consider two frequencies and show how the com- 
bination of information from both can significantly reduce the 
level of contamination. Clearly to make full use of the data set 
one should consider all the frequencies. Here we illustrate the 
method with just two but the generalization is trivial. 

6.1. Fisher matrix 

To illustrate how the cleaning works let us take a simple 
model for the data, 

aimiy) = af:"\v) + a{jy) + alTiy)- (38) 
The three terms are the 21 cm signal, foreground contami- 
nation, and detector noise. We will assume that we make 
Ni measurements at two separate frequencies, v\ and V2, so 
that the data vector is of the form x, = (a;,,, (;^i ),«;„, (j^2 ) ) with 
i=l,---Ni. We will assume that both the 2 1 cm signal and the 
detector noise are uncorrelated between the two frequencies 
while the foregrounds have a correlation coefficient /, very 
close to unity. In that case the correlation matrix of the data 
is 



(36) 



21c 



10 
01 



10^ 
01 



(39) 
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where (3 characterizes the frequency dependence of the fore- 
grounds and Cf gives the power spectrum of the noise, which 
for simplicity we assumed equal in both channels and uncor- 
related between the different Ni measurements (see §7). 

We will assume that the fluctuations are Gaussian and that 
we wish to estimate the three parameters of the model, p = 
(c/,/3,C,^ '''"), simultaneously. We can calculate the expected 
error bars from the Fisher matrix JF, 



-Tr 
2 



I dC 



(40) 



where h, 12 run over the three parameters. The inverse of 
the Fisher matrix gives the expected covariance matrix of the 
recovered parameters (see Tegmark, Taylor & Heavens 1997 
for a summary of the Fisher matrix technique). In particular, 
the error in the recovered 21 cm power spectrum is simply 



■ -^3,3 



2_ 

Ni 



2/3 



(1+/3) 



(1-/)C/ 



2/?(l -/)C/ 



(41) 



where we have assumed Cj S> C}^™ ,Cf . 

We see that the foreground power spectrum is suppressed 
by a factor 1—7; so, as long as (1 — I)c{ < cf^'^'" the noise 
introduced by foregrounds can handled with this technique. 

6.2. Another derivation 

We can obtain the same result as above by considering the 
following problem. Assume we measure the a/,„'s at two dif- 
ferent frequencies and call the results where ; runs over 
the number of observations. The data is intrinsically of the 
form 



Xi = fi + ei, 



(42) 



where f, is the part coming from the correlated foreground 
contribution, e, has the 2 1 cm signal plus the noise contribu- 
tion at the first frequency and Sj has the 21 cm signal, noise 
and an additional contribution from the uncorrelated part of 
the foregrounds. Just by diagonalizing the covariance ma- 
trix above one can show that this uncorrelated component has 
variance 2/3(1 —I)cf, to lowest order in (1 — /). For conve- 
nience we have introduced /i, = Si — /3'/^e,-. 

The (x, . y,) data fall on a straight line with unknown slope 
(v^) that we need to determine. The 21 cm signal is encoded 
in the deviations of the data from a perfect line. We can de- 
termine the slope by writing a simple for the best fit line. 



(43) 



where b is the slope that we are trying to determine. We can 
easily minimize x~ with respect to b. The information about 
the 21 cm fiuctuations is encoded in the value of at this 
minimum, for which we obtain: 



(44) 



Note that in terms of the underlying model variables. 



(45) 



We can take averages over the uncorrelated component /i, 
(which in a sense is the only random variable in this approach) 
to obtain 

(X'(^'™«)) = (A^-1)^^- (46) 
Thus we could use 

^=iV~r^'(^™«) (47) 

as the estimator for the variance, which contains the 21 cm 
signal. The mean and variance of this estimator are 



■■{1+13) 
N-l 



2/3 



2/3 



(1-/)C/ 



(l-/)C/(4&) 



This is almost the same as we obtained earlier The differ- 
ence can be traced to the fact that to estimate Q^'™ from S we 

need to subtract the contribution proportional to (1 — I)c{ ; 
this adds an extra piece to the variance that is second order 
in (1 — /)C/ . The important point, however, is that the fore- 
ground term appears only in the form (1 — I)cf. 

6.3. The cleaned foreground signal 

We can now illustrate how the technique we propose re- 
duces the importance of contamination from unresolved point 
sources. Figure 5 again shows the power spectrum at z = 10, 
with and without the extra power from reionization (the other 
curves in this figure are described in §7). We also show 
(1 - I)Cj for the "intermediate" point source model of Di 
Matteo et al. (2002) (dot-dashed line), assuming that point 
sources with 5 > 0.1 mJy have been removed. For the corre- 
lation coefficient of the maps produced by the point sources 
we took 1 — / = 8 X 10^^. This corresponds to the following 
choices: = 0.3, ln(i/2/t^i) = 1.3 x lO'^, {SC/(Jc? = 0.1, 
7 = and dlnCi/dln( = 0. This choice of 5( is consistent 
with the results of Cohen et al. (2003). Note that we have 
assumed the Poisson term to be negligible, as argued by Di 
Matteo et al. (2002), and we have neglected variations in the 
clustering length with spectral index, although we do include 
imperfect correlations between sources with different (. As 
the figure shows, the expectation is that once the frequency 
information is used, the contamination becomes significantly 
smaller than the signal we wish to measure. 

Oh & Mack (2003) noted another problem related to simple 
contamination: if the beamsize changes with frequency, the 
number of foreground sources in a given beam also changes 
with frequency. In our language, an interferometer baseline 
samples shghtly different /-modes at different frequencies. 
The most obvious solution is either careful beamsize control 
with frequency or fine coverage in Fourier-space coverage, 
both of which must be determined during the interferometer 
design. If this is not possible, an alternative method is to note 
that this "leakage" contamination is also highly correlated in 
frequency- space and can be removed with techniques similar 
to ours (see also Gnedin & Shaver 2003). 

Note that in our estimate of the foreground removal we have 
assumed that the 21 cm signal between the two frequency 
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Fig. 5. — Observability of the estimated 21 cm fluctuation signal. The solid 
black and dotted green curves in each panel are the angular power spectrum 
of the signal at z = 10 (with Ai/ = 0.2 MHz) for our toy model and for a 
fully neutral medium. The long dash-dotted curve in the top panel shows 
the estimated point source foreground signal (see text). The other curves 
show sensitivity estimates for three proposed experiments: PAST (blue short- 
dashed), LOFAR (red long-dashed), and SKA (magenta dot-dashed). The 
top panel shows the noise power spectrum while the bottom panel shows the 
eiTor in the estimated power spectrum. We assume four weeks of continuous 
observations for each experiment. 



channels is completely uncorrelated [i.e. the covariance ma- 
trix of the 21 cm signal is diagonal in equation (39)]. In 
reahty, the off-diagonal elements of the matrix will contain 
the cross-correlation of the two frequency bins, Cl^™(y\ , vj). 
The recovered signal is proportional to [1 ~ Cf'™\v\^V2)\. 
This essentially provides a minimum bandwidth for fore- 
ground removal using this technique: we need the radial sep- 
aration of the two frequency bins to exceed the typical cor- 
relation length. In the example we show, we have chosen 
L(Ai^) > R, so the true signal would decrease by a small 
amount. However, the method we have presented uses only 
two neighboring frequency channels to remove the contami- 
nants. In reality, the full bandwidth of the observation should 
be used in estimating the smooth component of the spectrum 
at each point on the sky, which will considerably improve the 
removal algorithm. This will allow the single-frequency sig- 
nal to be isolated without significant loss from correlations 
between neighboring frequencies. We conclude that point- 
source foregrounds do not present a significant problem for 
these measurements provided that they are smooth in fre- 
quency space. 



7. DETECTABILITY 

In this section we will examine the prospects for detecting 
the 21 cm fluctuations. We will follow the notation of White, 
Carlstrom, Dragovan, & Holzapfel (1999) where the formal- 
ism used to analyze data for CMB interferometers such as 
the Degree Angular Scale Interferometer (DASl) and the Cos- 
mic Background Interferometer (CBl) was presented. We note 
that similar sensitivity estimates have been made for the Gi- 
ant Metrewave Radio Telescope by Bharadwaj 8c Sethi (2001) 
and Bharadwaj & Pandey (2003). 



The measured flux in a visibility is 
dT 



y(u) 



d^n ^7i(n)A(n) e 



27r/u-n 



(49) 



where A(n) is the primary beam and dB^/dT converts tem- 
perature to flux. In the Rayleigh-Jeans part of the spectrum, 
dB^/dT = Iks/y?- The Fourier wavenumber u is related to 
Z by M = 1/2tt. We can express (57], (n) and A(n) in terms of 
their Fourier components. 



A(n): 



d^u5Tb{-a)e-^'""" 
d^uA{vi)e-^'''"-^. 



(50) 



The Fourier components 5Th{vL) are essentially the same as 
the fl/m's introduced earlier except that the spherical harmonic 
decomposition has been replaced by a Fourier decomposition, 
valid only over small patches of the sky that can be taken to be 
flat. The variance of 5Th{vL) is given by the power spectrum, 

{6Th{\ix)6Tb{n2)) = (5^(ui +U2)C/=2^„,. (51) 

In terms of these Fourier variables the variance of the temper- 
ature is 



-'I—Ittu 



dl- 



(52) 



which can be compared to the exact formula, 

- (2/+l)C, 



I 



47r 



(53) 



We now calculate the averaged value of the square of the 
observed visibilities in terms of the power spectrum, where 
the average is over an ensemble of possible skies, 

(|V(u)|2)= (^^) Y |A(U-U')PQ.2.„' 



.'m2 



(54) 



If the visibility is observed for a time f,, the averaged noise 
squared in each visibility is given by (Rohlfs & Wilson 2000) 

2 



/ 2kBT,y, 



1 



Ai/f,. 



(55) 



where Ti,, is the system temperature, Ai' is the bandwidth, 
and Aciish is the area of each individual antenna in the array. 

We can compare equations (54) and (55) to define the power 
spectrum of the noise. 



1 



(56) 



MshdBjdT J Avt,jd^u' |A(u-u')p 

The result of the integral in the denominator depends on the 
shape of the primary beam (i.e. the beam of the individual 
dishes). To get an approximate answer we can use the fact that 
A(u) is different from zero in an area d^u and has to integrate 
to one, so J d^u' |A(u — u')p ~ l/d^u. Moreover, the size of 
the primary beam and thus d^u is directly related to the area 
of the dishes. We can approximately use A^/,,/, = X^d^u. The 
power spectrum of the noise then becomes 



Tl (27r)2 



Avtvd'^u Ah'tyd^l 



(57) 
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This is equation (17) of White et al. (1999) with the mapping 
I = liTu. Perhaps an easier way to understand this equation 
is to note that after a time ty the noise in the Fourier space 
pixel corresponding to the observed visibility is simply = 
r^j/Ai^fy. The noise can be expressed in terms of a power 

spectrum using = d^uCi = d^l Ci / (27r)^. 

Interferometric CMB experiments such as DASI and CBI 
were conducted at much higher frequencies than those in 
which we are interested, so the arrays were small enough that 
they could be rotated to compensate for the Earth's rotation. 
Consequently, each pair of antennae could integrate for an ar- 
bitrarily long time on a single Fourier component of the sky. 
This is not feasible in the present case because the distances 
involved are much larger: the arrays must have baselines on 
the order of a kilometer. We therefore need to calculate the 
fraction of the total observing time to that any given baseline 
is being observed. This fraction of time will not be uniform 
across the Fourier plane, and the details will depend on the 
element configuration. 

We will make a simple estimate here assuming that the 
Fourier coverage is roughly uniform. For a particular max- 
imum separation of the antennae the interferometer will cover 
Fourier space up to a maximum Imax- Thus, owing to the 
Earth's rotation, a region of area tt/^^^^ of Fourier space will 
be covered. At any given instant, however, only an area 
Npairsd^l is being observed. Thus each visibility will be ob- 
served roughly for a time 
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where Npairs = Ndish {Ndish — 1 ) /2 and where we have assumed 
Ndish ^ 1- Combining equations (57) and (58) we obtain 
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We can think of the array as a big telescope with diameter 
D, large enough to make a measurement of mode Imax which 
covers a total area A,o,a/. However, only a fraction of that area 
is covered with telescopes, NdishA-dish- That covering fraction 
can also be expressed in terms of Imax and <fl. 
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In terms of fcover the noise power spectrum is. 
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If one is interested in achieving maximum sensitivity at a 
particular scale I with a fixed number of elements of a given 
size, it is best to pack the elements as close as possible be- 
cause Imax oc but f cover oc . This is achieved when the 
/ of interest is close to Z„,a^. The signal from the "bubbles" 
is located somewhere in the range / ~ 1000— 10000 so one 
needs an array of size D = /A/27r ~ 300m - 3 km to observe 
it. 

For realistic arrays, the distribution of telescopes will not 
be uniform, so the coverage in Fourier space will vary with I. 
For example, the array may have a core at the center where 
telescopes are closely packed and a more dilute configuration 



at larger separations. Thus the covering fraction for the /'s 
measured by the core of the array will be much larger than for 
the higher /'s. We can introduce a function /(/) that encodes 
the geometry of the array, in terms of which 
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For our simple case of uniform Fourier coverage, 

f(i) = f la^si 
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The system temperature at these frequencies will be domi- 
nated by the sky brightness temperature so the figure of merit 
to compare different experiments is simply /(/). We take T^ys 
to be roughly 200 K, so the noise power spectrum is approxi- 
mately 
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Thus, if we want an experiment to make a map with good sig- 
nal to noise in a matter of weeks it needs to have /(/) ~ 0.1 
on the scales of interest. Note that this noise estimate does not 
take cosmic variance into account. Because each observation 
samples only a finite region and must ultimately be compared 
to a statistical model of the Universe, at a certain point (when 
the signal to noise is of order unity in each Fourier mode), 
the experimental power becomes Umited by the finite field of 
view and one is better off increasing the area of the sky being 
covered rather than going deeper on the same spot. We have 
quoted the answer for a bandwidth of order 0.4 MHz. There 
is little to be gained by making the bandwidth much larger be- 
cause by doing so one enters the Limber regime on the angular 
scales of interest (see Figure 4 and discussion thereof). Once 
in the Limber regime both signal and noise scale as 1 / Av. 
The exact crossover into the Limber regime will depend on 
the sizes of the bubbles. On the other hand, if the bandwidth 
is much smaller than the typical correlation length of the 21 
cm features, the signal will be difficult to detect as it will be 
confounded by the foregrounds. 

If one is interested in a statistical detection of the power 
spectrum rather than imaging, the required observing time is 
significantly reduced simply because one can make several 
estimates of the power spectrum on scales smaller than the 
total field of view. The error in the power spectrum estimate 
is [equation (41)] 



^2 1 cm 



(65) 



The signal to noise increases by a factor, yjNi/2 ^ 
V Ndish / fcoZ ' r ~ V^totai/Adish- In terms of the /-modes sam- 
pled, ^JNiJt. ~ Ijlmin, where corresponds to the total an- 
gular size of the field of view. 

Moreover, we have included only one frequency channel 
in our estimate (after foreground subtraction). If the signal 
that one is seeking has a frequency width 1 MHz, stacking 
channels adds even more statistical power: the number of es- 
timates of goes like the number of channels used in the 
estimate. For example. Figure 3 shows that the power spec- 
trum of our toy model changes relatively little over Az 1, 
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corresponding to a frequency width of ~ 10 MHz. Thus, 
about 25 channels could be stacked without losing too much 
redshift information. 

Several experiments are now being designed to have the 
capabiUty to measure the 21 cm signal. One is a pro- 
posed dedicated experiment called the Primeval Structure 
Telescope^ (PAST). This instrument will have an effec- 
tive area Ndish^dish ~ lO'^m^ concentrated in a diameter 
D 1^ 2 km (l^ 5000). For that configuration /(/) ~ 
NdisiAdish/D^ (Imax/l) ~ 0.0024 l„iax/l- Thus thc instrument 
would require long integrations or averages over many fre- 
quency channels to detect the expected statistical signal. Note 
that some long basehnes would also be needed to be able to 
remove point source contamination. 

Detecting cosmic 21 cm emission is also one of the major 
science goals of the Low Frequency Array^ (LOFAR) and the 
Square Kilometer Array'" (SKA). LOFAR will have a total 
effective area of about 2 x 10^ m^ with approximately 25% of 
that area concentrated in a compact array of D ^ 2 km. For 
this core, Imax ~ 5000 and /(/) ~ 0.016 lm<a/l- The design of 
the SKA has not yet been fixed. Current plans call for ~ 20% 
of the array elements to lie in a core oi D ^ \ km and ^ 50% 
to he within a region of Z) ~ 6 km. For the inner region, Imca ~ 
2500 and /(/) ~ 0.25 Imwc/l and for the outer one Imax ^ lO'' 
and f{T) ^ 0.018 Imax/l- Both of these instruments also have 
the advantage of very long baselines (hundreds or thousands 
of kilometers) that will help with point source removal and 
control of systematics. 

Figure 5 shows some estimated sensitivity curves for each 
of these instruments. To construct the curves, we assumed 
one month of continuous observing on a single field of view. 
We assumed a 100 square degree field of view for each of the 
experiments. We have also made some assumptions about the 
antenna distributions in these experiments. For LOFAR, we 
took (25%, 50%) to have baselines smaller than (2, 12 km). 
For SKA, we took (20%, 50%, 55%) within (1,6, 12 km). We 
assumed that /(/) varies smoothly between these points and 
uniform Fourier space coverage within the core. Note that our 
sensitivity estimates are only approximate and depend on the 
Fourier space coverage of the array as well as the correlation 
procedure. 

The top panel shows the noise power spectrum; compar- 
ison to the power spectrum of the 21 cm fluctuations gives 
the signal-to-noise value for each measured visibility. This 
panel therefore shows the appropriate sensitivity for making a 
map. Clearly, creating images with high signal to noise on 
arcminute scales will be difficult and require large collect- 
ing areas (on the order of a square kilometer). Note that the 
slope of the sensitivity curves depend on the antenna config- 
uration; for uniform Fourier coverage, ST„ oc I. Configura- 
tions in which the coverage increases at smaller separations 
(i.e. the array's covering fraction increases toward the center) 
have steeper slopes. 

Fortunately, as noted above, a statistical measurement of the 
power spectrum is considerably easier. In the bottom panel 
we show the error on the estimated power spectrum in loga- 
rithmic Z-bins: this is simply the noise power spectrum multi- 

phed by A^; ' ~ Imm/l- The large field of view of SKA allow 
rather precise measurements of the power spectrum on scales 
from < 1 arcmin to 1°. Note also that we show the errors 

* See http://astrophysics.phys.cmu.edu/~jbp for details on PAST. 
' See http://www.Iofar.org for details on LOFAR. 

See http://www.skatelescope.org for details on the SKA. 



in the individual frequency channels. If many channels are 
stacked together, the errors will decrease by the square root of 
the number of channels (ignoring correlations between chan- 
nels). In this way, PAST and especially LOFAR could also 
make significant detections of the power spectrum 

8. CONCLUSIONS 

The recent successful efforts to map anisotropies in the 
cosmic microwave background have determined the global 
properties of the Universe to high precision (e.g. Spergel 
et al. 2003). When combined with the highly successful 
paradigm for the hierarchical growth of structure (e.g. White 
& Rees 1978), these results imply that the evolution of the 
Universe on large scales is now relatively well-understood. 
However, on the smaller scales characteristic of individual ob- 
jects, many uncertainties remain owing to our ignorance of the 
nature of dark matter, the lack of a full physical model for the 
origin of primordial density fluctuations, and our poor under- 
standing of galaxy formation. Determining the properties and 
consequences of the first luminous objects at z 15 — 30 may 
help to clarify these issues. 

The evolution of the ionized part of the IGM at high red- 
shifts provides a fossil record of the Universe at these times. 
In principle, the physical state of this diffuse gas constrains 
when and where the first luminous objects formed as well as 
the nature of the sources responsible for reionization. The rel- 
atively large electron scattering optical depth obtained from 
the WMAP observations is indicative of a complex reioniza- 
tion history but, by itself, does not provide unambiguous an- 
swers to the remaining questions. 

In this paper, we have argued that fluctuations in 21 cm 
emission from the IGM can be used to measure the three- 
dimensional distribution of neutral hydrogen at high redshifts. 
Measurements of the angular power spectrum at different fre- 
quencies can be used to mitigate contamination by foreground 
sources that would otherwise overwhelm the 21 cm signal. 
Our approach is similar to that employed in analyses of CMB 
anisotropies, but is more general because of the frequency- 
dependent nature of redshifted 21 cm emission, making it pos- 
sible to use 21 cm fluctuations study the evolution of reioniza- 
tion. We note here that the most basic kind of experiment is 
to seek an "all-sky" signature corresponding to a global phase 
transition in the neutral hydrogen (which could be reioniza- 
tion, reheating, or the onset of Lya coupling; Shaver et al. 
1999). While valuable (especially because they have much 
less stringent sensitivity requirements), these measurements 
only provide the most basic information. Moreover, they are 
subject to severe foreground contamination (Gnedin & Shaver 
2003) because the large scale signal varies relatively smoothly 
with frequency. Thus frequency differencing will not be as ef- 
fective as with small-scale fluctuations. The small-scale fluc- 
tuations we have described will therefore ultimately be a bet- 
ter route to pursue. 

Using a simple conceptual model for the morphological 
evolution of ionized regions, we have demonstrated how 
reionization imprints characteristic features into the angular 
power spectrum of 21 cm fluctuations. Our formulation is 
general, but for illustrative purposes we have adopted sim- 
plifying assumptions to emphasize salient features of our 
methodology. For example, in our toy model of reionization, 
we have ignored redshift space distortions in the 21 cm sig- 
nal and neglected correlations between the density and neu- 
tral fraction fields. In future work, we will investigate the 
impact of these effects explicitly using semi-analytic methods 
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and numerical simulations to show how the formalism can be 
adapted to account for these complications. 

In principle, measurements of the type we propose can be 
used to distinguish between various evolutionary histories. 
Models with multiple epochs of reionization (e.g. Cen 2003; 
Wyithe & Loeb 2003) lead to a behavior in which the volume 
fraction of ionized gas in the universe shows a complex de- 
pendence on redshift (e.g. Figures 8 & 9 of Sokasian et al. 
2003b), unlike single episodes of reionization where the trend 
is simpler (e.g. Figure 5 of Sokasian et al. 2003a). These dif- 
ferences will be imprinted onto the frequency dependence of 
21 cm fluctuations of the IGM and can be discerned using the 
approach described here. In fact, the simple model described 
in §4 is only part of the story available to us through 2 1 cm 
observations. For example, the fluctuations constrain the ther- 
mal history of the IGM as well as the ionization history (see 
§2). If the earliest ionizing sources have hard spectra (such 
as quasars), we would expect the IGM to be heated rapidly, 
while if cool, low-mass stars are responsible for reionization, 
the heating would occur much closer to overlap. Another way 
to look at this is that the reionization history determines how 
much information is available to us through 21 cm fluctua- 
tions. 

For instance, in a scenario with rapid heating, we will have a 
long epoch where density variations dominate the signal. The 
power spectrum of the 21 cm fluctuations is a direct tracer 
of the matter power spectrum. Thus, detailed measurements 
of this signal could provide invaluable constraints on the pri- 
mordial spectrum on small scales which could tightly con- 
strain the physics of the early Universe. In particular, because 
one gets many independent maps by varying the frequency, 
the constraints on the power spectrum obtained in this way 
could be significantly better than those coming from the CMB 
primary anisotropics. Also, 21 cm fluctuations allow one to 
probe the power spectrum at z ~ 10-20, an epoch inacces- 
sible to both the CMB and galaxy surveys and a relatively 
large lever with which to study the growth of fluctuations. On 
the other hand, the 21 cm signal depends on complex physi- 
cal processes (see §2) and interpreting the results will require 
careful modeling. 

Differences in the evolutionary state of the ionized IGM 
also contain information about the matter power spectrum on 
small scales. For example, in cosmological models with re- 
duced small-scale power, the halos hosting star-forming re- 
gions form late, delaying the ionizing effects of the first lu- 
minous objects (e.g. Somerville et al. 2003; Yoshida et al. 
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2003b, c). Thus, universes with a large component of warm 
dark matter or those in which the matter power spectrum has 
a running spectral index should exhibit lower amphtude fluc- 
tuations in 21 cm emission from higher redshifts than in con- 
ventional ACDM models. 

A detailed study of 21 cm fluctuations can also constrain 
the nature of the sources responsible for reionization. In some 
scenarios, it is conjectured that reionization occurs "outside- 
in," affecting voids first and high density regions later (e.g. 
Miralda-Escude et al. 2000). This progression would be ex- 
pected if the sources are rare, but bright. Alternatively, reion- 
ization could proceed in the opposite sense, "inside-out," par- 
ticularly if the sources are numerous, but faint (e.g. Gnedin 
2000; Sokasian et al. 2003a). The former possibility would 
apply if quasars were the primary sources of ionizing radia- 
tion, while stars in low-mass galaxies would be more relevant 
to the "inside-out" scenario. The morphological appearance 
of the ionized gas is different in these two cases, and these 
differences would be reflected in the angular power spectrum 
of 21 cm fluctuations and how this quantity varies with fre- 
quency. In particular, the size of the H II regions at a given 
xh constrains the number density of ionizing sources. Also, 
if voids are ionized first the density and xh fields will be cor- 
related (i.e., the neutral fraction falls first in regions that are 
already underdense), while the opposite is true in an "inside- 
out" scenario. Thus the amplitude of the brightness temper- 
ature fluctuations contains information about the process of 
reionization. 

As we have argued in §7, the technological requirements 
for detecting 21 cm fluctuations from cosmic gas at high red- 
shifts, while demanding, are within reach. In the near future, 
it is likely that measurements of the power spectrum of 21 cm 
fluctuations will reveal the physical state of the Universe at 
an epoch that is currently inaccessible to other observational 
probes. 
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